# Set parent dir for visium experiments containing folders of spaceranger output with names
visium.dir <- "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RRRM2_ISST/Hippocampus"
# Path for where the plots will be exported to
export_path = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RRRM2_ISST"
samples <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'filtered_feature_bc_matrix.h5')
spotfiles <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_positions.csv')
imgs <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_hires_image.png')
json <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'scalefactors_json.json')
section.name <- samples
section.name <- gsub(paste0(visium.dir, "/"),"", gsub("/filtered_feature_bc_matrix.h5","",section.name))
infoTable <- data.frame(section.name, samples, spotfiles, imgs, json, stringsAsFactors = FALSE)
VlnPlot(se_brain, features = "nFeature_RNA", ncol = 1, group.by = "section.name", pt.size = 0)
ST.FeaturePlot(se_brain, features = "nFeature_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)
VlnPlot(se_brain, features = "nCount_RNA", ncol = 1, group.by = "section.name", pt.size = 0)
ST.FeaturePlot(se_brain, features = "nCount_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)
se_brain$percent.mt <- PercentageFeatureSet(se_brain, pattern = "^mt")
se_brain$percent.ribo <- PercentageFeatureSet(se_brain, pattern = "^Rpl|^Rps")
VlnPlot(se_brain, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.ribo"), ncol = 2, pt.size = 0)
Run PCA here on the non-filtered data and then plot it
p1 <- DimPlot(se_brain_merged_notfilt, group.by = "section.name")+ ggtitle("Plot of the non-filtered data")
p2 <- DimPlot(se_brain_merged_notfilt)
p1|p2
coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
DimPlot(se_brain_merged_notfilt, cols = coldef)
enids <- read.table(file = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RR3_test/rr3_ref/mm10_genes.tsv", header = T, stringsAsFactors = T)
enids <- data.frame(apply(enids, 2, as.character), stringsAsFactors = F)
rownames(enids) <- enids$gene_id
se_brain <- SubsetSTData(se_brain, expression = nFeature_RNA > 200)
staffli <- se_filtered@tools$Staffli
se_brain_merged <- merge(se_section_list[[1]], y = c(se_section_list[2:length(se_section_list)]), merge.data = T)
VariableFeatures(se_brain_merged) <- features
Plot data before and after harmony
se_brain_merged <- RunPCA(se_brain_merged)
## PC_ 1
## Positive: Nrgn, Camk2n1, Slc17a7, Olfm1, Nptxr, Snap25, Cck, Chn1, Ppp3ca, Hpca
## Ctxn1, Calm2, Atp1a1, Rtn1, Snca, Calm1, Lamp5, Camk2a, Enc1, Ppp3r1
## Gpm6a, Ncdn, Mef2c, 1110008P14Rik, Aldoa, Hpcal4, Ptk2b, Cpne6, Stx1a, Slc30a3
## Negative: Plp1, Mbp, Mobp, Ptgds, Apod, Cnp, Trf, Mal, Fth1, Mag
## Cryab, Qdpr, Plekhb1, Scd2, Cldn11, Car2, Mog, Sept4, Csrp1, Gatm
## Dbi, Apoe, Sparc, Ndrg1, Ppp1r14a, Dbndd2, Ermn, Bcas1, Pllp, Mt1
## PC_ 2
## Positive: Pcp4, Prkcd, Ramp3, Tcf7l2, Tnnt1, Rora, Adarb1, Ntng1, Plekhg1, Uchl1
## Cplx1, Pdp1, Ptpn4, Ndrg4, Slc17a6, Ccdc136, Plcb4, Zic1, Ptpn3, Amotl1
## Shox2, Rgs16, Nsmf, Pcp4l1, Patj, Kcnc2, Synpo2, Atp2a2, Nrip3, Gabra4
## Negative: Nrgn, Fth1, Ttr, Plp1, Cnp, Mal, Tmsb4x, Mobp, Enpp2, Apod
## Ddn, Clu, Ctxn1, Cryab, Trf, Gfap, Mbp, Aplp1, Cldn11, Hpcal4
## Cst3, Hpca, Camk2a, Penk, Mag, Ptgds, Camk2n1, Cplx2, Psd, Nptxr
## PC_ 3
## Positive: Nrgn, Plp1, Mbp, Fth1, Snap25, Cck, Mobp, Slc17a7, Camk2n1, Pcp4
## Prkcd, Mal, Scn1b, Cnp, Mag, Cryab, Cplx1, Slc24a2, Trf, Hpca
## Chn1, Sept4, Rasgrp1, Ppp3ca, Adcy1, Qdpr, Car2, Dbndd2, 1110008P14Rik, Mog
## Negative: Nnat, Resp18, Hap1, Nap1l5, Ndn, Zcchc12, Calb2, Gaa, Gpx3, Sparc
## Tmem130, Ahi1, Baiap3, Ly6h, Wdr6, Gap43, Cartpt, AW551984, Scg2, Nrsn2
## Agt, Fxyd6, Gprasp2, Pnck, Pgrmc1, Vat1l, Tac1, Ngb, Bex1, Dlk1
## PC_ 4
## Positive: Plp1, Mbp, Fth1, Resp18, Mobp, Nap1l5, Cnp, Stmn1, Mag, Mal
## Qdpr, Ndn, Stmn3, Sept4, Tubb4a, Cryab, Tmsb10, Gap43, Uchl1, Mog
## Tmem130, Trf, Ubb, Zwint, Dbndd2, Bex2, Zcchc12, Hap1, Gatm, Gapdh
## Negative: Apoe, Cst3, Clu, Slc1a2, Hbb-bs, Camk2a, Mt1, Atp1a2, Igfbp2, Glul
## Igf2, Hba-a1, Hba-a2, Ttr, Ddn, Mt2, Pltp, Ptgds, Lars2, Id3
## Mgp, Vtn, Rbp1, Mt3, Fxyd1, Hbb-bt, Vim, Gfap, Dbi, Ifitm3
## PC_ 5
## Positive: Penk, Pcp4, Ncdn, Ppp1r1b, Gpr88, Pde1b, Hpca, Rgs9, Wipf3, Wfs1
## Gng7, Tac1, Ddn, Pde10a, Fkbp1a, Prkcd, Adora2a, Adcy5, Necab2, Ppp3ca
## Ptk2b, Prkcg, Camkv, Tmsb4x, Cpne6, 2010300C02Rik, Lrrc10b, Cnih2, Crym, Atp2b1
## Negative: Snap25, Pvalb, Camk2n1, Lamp5, 1110008P14Rik, Nrgn, Atp1a1, Stmn1, Kcnab3, Vsnl1
## Cabp1, Mef2c, Scn1b, Osbpl1a, Stx1a, Ttr, Serpini1, Cplx1, Nap1l5, Syt2
## Vamp1, Egr1, Igfbp6, Olfm2, Hbb-bs, Gabra1, Gnas, Sncb, Arc, Nrsn1
se_brain_merged <- RunUMAP(se_brain_merged, dims = 1:35, reduction = "pca")
## 09:49:38 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:49:38 Read 41301 rows and found 35 numeric columns
## 09:49:38 Using Annoy for neighbor search, n_neighbors = 30
## 09:49:38 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:49:41 Writing NN index file to temp file /var/folders/t8/1flxfkx50cb3ycdyf5y0gb0h0000gn/T//RtmpDtghuS/file12ee3ad9ebd8
## 09:49:41 Searching Annoy index using 1 thread, search_k = 3000
## 09:49:50 Annoy recall = 100%
## 09:49:50 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:49:52 Initializing from normalized Laplacian + noise (using irlba)
## 09:49:53 Commencing optimization for 200 epochs, with 1892064 positive edges
## 09:50:18 Optimization finished
p3 <- DimPlot(se_brain_merged, group.by = "section.name")+ ggtitle("Plot of the filtered data before integration with Harmony")
p4 <- DimPlot(se_brain_merged)
p3|p4
coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
DimPlot(se_brain_merged, cols = coldef)
p7 <- DimPlot(se_merged_harmony, group.by = "section.name")+ ggtitle("Plot of the filtered data after integration with Harmony")
p8 <- DimPlot(se_merged_harmony)
p7|p8
coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
DimPlot(se_merged_harmony, cols = coldef)
se_merged_harmony@tools$Staffli <- staffli
ST.FeaturePlot(se_merged_harmony, features = "seurat_clusters", dark.theme = F, pt.size = 2.3, show.sb = F, ncol = 4, cols = coldef)